Spin Flips and Precession in Black-Hole-Binary Mergers 
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We use the 'moving puncture' approach to perform fully non-linear evolutions of spinning quasi- 
circular black-hole binaries with individual spins unaligned with the orbital angular momentum. 
We evolve configurations with the individual spins (parallel and equal in magnitude) pointing in 
the orbital plane and 45° above the orbital plane. We introduce a technique to measure the spin 
direction and track the precession of the spin during the merger, as well as measure the spin flip in 
the remnant horizon. The former configuration completes 1.75 orbits before merging, with the spin 
precessing by 98° and the final remnant horizon spin flipped by ~ 72° with respect to the component 
spins. The latter configuration completes 2.25 orbits, with the spins precessing by 151° and the final 
remnant horizon spin flipped ~ 34° with respect to the component spins. These simulations show 
for the first time how the spins are reoriented during the final stage of black-hole-binary mergers 
verifying the hypothesis of the spin-flip phenomenon. We also compute the track of the holes before 
merger and observe a precession of the orbital plane with frequency similar to the orbital frequency 
and amplitude increasing with time. 

PACS numbers: 04.25.Dm, 04.25. Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 



There is widespread interest in understanding the dy- 
namics of the last orbital stages of generic black- hole bi- 
naries (i.e. binaries with randomly aligned spins and un- 
equal masses). These last few orbits before the binaries 
merge involve highly-nonlinear gravitational interactions 
leading not only to large amounts of gravitational radi- 
ation leaving the systems, but also to intricate coupling 
effects, particularly those involving the spins of the two 
component black holes. Questions of how the black holes 
orient their spins with respect to the (instantaneous) or- 
bital plane and the extent to which the black holes spin- 
up during the last orbital stages are of great astrophysical 
interest. Spinning black holes are believed to be the en- 
gines of active galactic nuclei and quasars; the efficiency 
with which infalling matter is converted into radiation de- 
pends on the energy of the innermost stable orbit and is 
greatest for rapidly-spinning holes [l], Q • The energetic 
jets observed in many AGN and stellar-mass accreting 
systems are believed to be launched perpendicularly to 
the inner accretion disk, hence parallel to the spin axis 
of the accreting hole Changes in the spin direction 
are therefore potentially observable. A number of active 
galaxies exhibit semi-periodic deviations of the jet direc- 
tions from a straight line [H, suggestive of precession of 
the accretion disk around the jet-emitting hole or geode- 
tic precession of the larger hole, either of which might be 
driven by torques from a second orbiting black hole [j| @ . 
About 15 radio galaxies show jets with apparently abrupt 



changes in jet direction, forming X-shaped patterns 
Jets of Seyfert galaxies often misalign from subkiloparsec 
to kiloparsec scales [|| . Some models [T3, E| attribute 
these spatial variations to a sudden re-orientation of the 
spin axis of the larger hole as it accretes a smaller hole. 
Fossil evidence that such re-orientations were common in 
the past is observed in the nearly random orientations of 
jets in disk galaxies with respect to the disk plane fl2l.[l3|. 

In the generic case, inspiral of a black-hole binary 
should induce a precession in the direction of the spin 
axis of either hole [3, and unless the mass ratio is ex- 
treme, the final spin orientation is dominated by the or- 
bital angular momentum of the binary, implying substan- 
tial re-orientation as long as the spin and orbital angular 
momenta are initially mis-aligned [lOj]. We refer to this 
jump in the spin direction of the remnant, with respect 
to the individual spins, as a spin-flip. Note that this 
definition of spin-flip does not require non-linear interac- 
tions; it is a simple consequence of the remnant acquiring 
most of the total angular momentum of the original bi- 
nary. Also note that the 'flip' is actually the difference 
in spins between two distinct types of objects: individual 
horizons and common horizons. 

Aside from their astrophysical interest, black- hole bi- 
naries are one of the primary targets of the earth-based 
gravitational wave observatories, such as LIGO [Hj], 
VIRGO 0, GEO600 [H, etc. These detectors are 
operating at unprecedented sensitivities. In particular, 
the LIGO interferometers are currently taking data at 
their design sensitivities. The initial LIGO interferome- 
ters might be able to observe two 20M Q inspiraling black 
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holes out to a distance of more than lOOMpc, and ad- 
vanced LIGO could see the same event up to a cosmolog- 
ical redshift of z « 0.4 [l8(. Predicted black-hole-binary 
gravitational waveforms will not only be of great help 
for the detection of this radiation using matched filtering 
techniques, but will also be essential for the interpreta- 
tion of the signals, determination of the event rates, and 
extraction of astrophysical parameters. This will be es- 
pecially important for the next generation of ground and 
space-based detectors, such as LISA [l9|,[2(| which should 
observe gravitational wave bursts from the mergers of su- 
permassive black holes in the centers of galaxies out to 
very high redshift s. 

Simulations of the last orbital stages of black-hole bi- 
naries require solving the fully-nonlinear General Rel- 
ativity field equations numerically on supercomputers. 
However, solving these non-linear equations proved to 
be quite difficult and the problem remained unsolved for 
over thirty years. However, last year two independent 
techniques were developed that broke through the bar- 
rier of the numerical instabilities to produce spectacu- 
lar results. In 2005 these two approaches were used to 
generate the gravitational waveforms from the last or- 
bit of non-spinning, equal-mass black-hole binaries. The 
first technique, which was developed by Pretorius [2l| . 
used a second-order formulation of the General Rela- 
tivity field equations in a generalized harmonic gauge 
(GHG), along with singularity excision in the interior 
of the horizons, adaptive mesh refinement (AMR) on a 
compactified space, and the addition of constraint damp- 
ing terms to the evolution equation. The second success- 
ful approach was developed a few months later by our 
group at UTB [22j and independently by the Numeri- 
cal Relativity group at NASA/Goddard (23|. This latter 
method uses a mixed flrst-and-second-order formulation 
of the General Relativity field equations known as BSSN 
system [24l l25l l2~6ll . in combination with the puncture 
formalism [23, [28| (without the need for singularity ex- 
cision). This technique differs from previous work with 
punctures in that the punctures are not fixed on the grid. 
In both versions of this 'moving puncture' approach mod- 
ifications to the standard I+log lapse and Gamma-driver 
shift conditions were introduced [H, US . 

Most of the groups in numerical relativity have now 
implemented one of the two approaches; the 'moving 
puncture' technique being the more popular. Given its 
technical simplicity and flexibility, this latter approach 
has been used to produce several interesting results rang- 
ing from the original simulations of non-spinning equal- 
mass binaries [3fl l3ll . l32l |33| to spinning equal-mass 
binaries [34l, 13511 and non-spinning unequal-mass bina- 
ries [36l l37l . |38| |. Notably, the moving puncture tech- 
nique has recently also been successfully implemented 
in the case of neutron-star — black-hole-binaries [3!| H3] , 
thus extending the impact of the technique on the nu- 
merical/ astrophysical relativity community. Noteworthy, 
the GHG approach has also been successfully used to 
study eccentric orbits in black-hole binaries with equal- 



mass and small corotation spins [3I[ . A generalized har- 
monic form of Einstein's equations has also been suc- 
cessfully used together with a dual coordinates method 
to evolve black-hole-binary spacetimes for several orbits 
prior merger [4l| . 

Alongside this impressive progress in computational 
and experimental relativity, the last few years have also 
seen significant mathematical progress in our under- 
standing of black holes in full non-linear General Relativ- 
ity. There is now a better understanding of the geometry 
and dynamics of trapped and marginally trapped surfaces 
using the quasi- local notions of trapping [42j , isolated [43[ 
and dynamical horizons [44] . Isolated horizons describe 
black holes in equilibrium in non-stationary spacetimes, 
and trapping and dynamical horizons describe the gen- 
eral time-dependent case. The applications of these ideas 
to classical and quantum black hole physics are too nu- 
merous to be described here, and we refer to [HI, H(| [47| 
for reviews and a more complete set of references. For 
our present purposes, we are mostly interested in cal- 
culating the angular momentum of an (approximately) 
axisymmetric horizon. The calculation of angular mo- 
mentum for isolated horizons is carried out using Hamil- 
tonian methods as described in [48|, H^, [5(| ■ The analo- 
gous Hamiltonian calculation for non-stationary trapping 
and dynamical horizons is given in 51]. Conservation 
and balance laws describing how the horizon mass and 
angular momentum change in response to infalling mat- 
ter / radiation are found in [5^, [Hf, HU, [5f| . The calcula- 
tion of the magnitude of black hole spin angular momen- 
tum tailored to numerical relativity is presented in [56j 
and more recently, [57| considers higher mass and angu- 
lar momentum multipole moments. In this paper we use 
this formalism primarily to compute the direction and 
magnitude of the spin angular momentum vector. 

This paper is organized as follows. In Sec.[TT]we discuss 
the Post-Newtonian predictions for spin and orbital plane 
precession. In Sec. IIIII we describe the techniques used 
to evolve the binary and measure the horizon spins. In 
Sec. IIVI we describe the initial data parameters for the 
binary configurations mentioned in the remainder of the 
text. In Sec. [V] we give a detailed description of the new 
results regarding spin and orbit precession. In Sec. lVII we 
discuss some of the implications of our results. Finally, 
in Appendix [X] we review past results from aligned-spin 
and non-spinning binaries. 



II. POST-NEWTONIAN ANALYSIS 

If the spins of the two black holes in a binary are not 
aligned with the total angular momentum, then the spin 
and orbital angular momenta will precess about the total 
angular momentum. Precession of the spin of the holes 
is produced by spin-orbit coupling and the spin-spin cou- 
pling. This effect has been studied in several papers by 
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means of the Post- Newtonian expansion [58|, [5 
dSi If , =» =? . ,„ 3 m,- 



(Ljv x Si)(2 



2 m. 



*5j" x Si 



+ 3(n • S^n x Si 



(1) 



where L^v = [a(x X v) is the Newtonian orbital angular 
momentum, x = X\ — X2, r = |a?|, v = dx/dt, n = x/r, 
ni = m-y + m 2 , [i = mim^/m, rj = fi/m, S = S± + S2, an 
over-dot denotes d/dt, and the j subscript denotes the 
other hole. Since the time derivative of the individual 
spin of the holes has the form ^ = !1, x Si this implies 
that the magnitude of each individual spin is conserved 
and only its direction will continuously change in time 
and rotate at a precession frequency given by 
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Sj + 3(n • Sj)n 



(2) 



Both the spin and orbital planes precess. This is due 
to the fact that, if the radiated angular momentum is 
neglected, the total angular momentum J = L + S is 

conserved. Hence the L = —S. Thus the orbital plane 
will precess at the same frequency as the total spin. The 
gravitational radiation reaction effect will generate a net 
loss of J that will actually produce an increase in the 
amplitude of the spin and orbital oscillations since the 
black holes will get closer; magnifying the spin-orbit cou- 
pling. In this paper, we go beyond the Post-Newtonian 
expansion and study these effects using full numerical 
evolutions. 

Precession occurs in the plane perpendicular to the to- 
tal angular momentum. We thus split the spin vector of 
each hole into components parallel to, and perpendicular 
to, the total angular momentum, i.e. S = S\\ + S±, where 

S\\ = {S- J) J and S± = S- (S- J) J. Note that the direc- 
tion (but not magnitude) of the total angular momentum 
does not change significantly between the start and end 
of the simulations. We define the total angle of preces- 
sion 9 p as the angle by which S± is rotated between the 
start and the end of the simulation. The total precession 
angle is then given by 



cos0 p = 



Sm ■ Si — (J ■ Sm)(J ■ Si) 

[s 2 M - (J ■ s M ) 2 ][s 2 - (J • § T y 



(3) 



where Si is the initial spin of one of the individual hori- 
zons, Sm is the spin of that individual horizon at the 
merger time, and J is the initial total angular momen- 
tum of the system. 

The relevance of processing spinning-black-hole bina- 
ries to data analysis has been stressed in several papers 
using the post-Newtonian templates. 'Spiky' templates 
have been considered to detect moderate massive galactic 
binaries in [Ml, H3| • Detection and post-Newtonian 



dynamics in processin g b inaries have also been exten- 
sively discussed in [63l 164 l65l l66l l67l . l68j , and the rele- 
vance to LISA observations has been discussed in 691. 



III. TECHNIQUES 

We use the Brandt-Briigmann puncture ap- 
proach [27| along with the TwoPunctures [70| 
and BAM_Elliptic [71( thorns to compute initial 
data. In this approach the 3-metric on the initial 
slice has the form 7^ = (ipBL + u) 4 3ab, where ipBL 
is the Brill-Lindquist conformal factor, Sab is the 
Euclidean metric, and u is (at least) C 2 on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
ipBL = 1 + Yn=i m i/(2 r i)i where n is the total number 
of 'punctures', rrii is the mass parameter of puncture i 
(rrii is not the horizon mass associated with puncture 
i), and n is the coordinate distance to puncture i. 
In all cases below, we evolve data containing only 



two punctures with equal puncture mass parameters, 
and we denote this puncture mass parameter by m p . 
We evolve these black-hole-binary data-sets using the 
LazEv I72I imp lementation of the moving puncture 
approach [22 . 23]. In our version of the moving puncture 
approach j22j we replace the BSSN 0, [H, [26| conformal 
exponent 0, which has logarithmic singularities at the 
punctures, with the initially C 4 field x = exp(— 40). 
This new variable, along with the other BSSN variables, 
will remain finite provided that one uses a suitable 
choice for the gauge. An alternative approach uses 
standard finite differencing of </> [23|. Note that both 
approaches have been used successfully by several other 

groups n, m m [zi izi iii . 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1+log lapse, a modified Gamma-driver shift 
condition f22l . |29( |. and an initial lapse a ~ ^1 The 
lapse and shift are evolved with (pt — (3 l di)a = —2aK, 
d t /3 a = B a , and d t B a = 3/4<9 t f Q - r\B a . These gauge 
conditions require careful treatment of x near the punc- 
ture in order for the system to remain stable [H, [30, lzl| • 
In Ref. [76| it was shown that this choice of gauge leads 
to a strongly hyperbolic evolution system provided that 
the shift does not become too large. For our version of 
the moving puncture approach, we find that the prod- 
uct aA^djip initially has to be C 4 on the puncture. In 
the spinning case, A lJ is 0(r 3 ) on the puncture, thus 
requiring that a oc r 3 to maintain differentiability. We 
therefore choose an initial lapse a(t = 0) = 2 / (1 + iI>bl) 
which is 0(r 4 ) and C 4 on the puncture and reproduces 
the isotropic Schwarzschild lapse at large distances from 
the horizons. The initial values of (3 l and B % are set to 
zero. 

Hannam et. al. 74j examine the smoothness of the 
evolved fields at late times at the puncture. They find 



that, in the case of Schwarzschild, x transitions from an 
initially C 4 field to a C 2 field at late times. Although we 
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require that the fields be initially C 4 , this late-time drop 
in smoothness does not appear to leak out of the horizon 
(which is consistent with the analysis in [74j]). 

We use a 'multiple transition' fisheye transforma- 
tion [3(| to push the boundaries to 200M, while main- 
taining a resolution of up to M/30 in the central region. 

We measure the magnitude S of the angular momen- 
tum of the horizons using our implementation of the al- 
gorithm detailed in [56]. The magnitude of the horizon 
spin is given by 

S = ±- / (<p a R b K ab )d 2 V (4) 

87r J AH 

where (p a is an approximate Killing vector on the hori- 
zon, K ab is the extrinsic curvature of the 3D-slice, d 2 V 
is the natural volume element intrinsic to the horizon, 
and R a is the outward pointing unit vector normal to 
the horizon on the 3D-slice; the sign of (p a is chosen so 
that S is positive. This algorithm for calculating S was 
initially meant to be applied to the case when the indi- 
vidual black holes are modeled as axisymmetric isolated 
horizons, which is valid when the two black holes are suf- 
ficiently far away from each other. The isolated horizon 
formalism is generalized to the dynamical case through 
the notion of a dynamical horizon [521 ] , and the formula 
for S remains valid under this generalization. 

Turning now to the direction of the spin angular mo- 
mentum vector, we first note that, in general, it seems dif- 
ficult to assign a unique coordinate independent 3- vector 
S to a spinning horizon. For example, we could take a 
normal Kerr spacetime and slice it non-axisymmctrically 
so that it becomes difficult to assign a spin 3-vector S 
to the black hole on these distorted 3D-slices. There is 
however a generalization which works. To see this, first 
note that every smooth cross section (with complete S 2 
topology) of a Kerr horizon is axisymmetric, no matter 
how distorted this cross-section is. This may seem some- 
what surprising at first glance, but it is a straightforward 
consequence of the fact that the null generators of the 
Kerr horizon have vanishing expansion, shear and twist; 
the axial symmetry vector projects to a symmetry of the 
2-geometry of the cross-section. Thus there exists a sym- 
metry vector ip a on this cross-section. The poles of the 
horizon are then defined to be the points where the axial 
symmetry vector ip a vanishes. From a spacetime perspec- 
tive, the locus of points on the Kerr horizon defined by 
ip a ip a — is a coordinate and gauge independent notion. 
These considerations remain valid on every axisymmetric 
isolated horizon. As long as we have a suitable axial vec- 
tor, we can similarly define the poles even for dynamical 
and trapping horizons. The poles exist whenever we can 
assign a (possibly approximate) axial symmetry vector 
ip a on the horizon. Of course, when the horizons become 
extremely distorted, it might happen that it is no longer 
approximately axisymmetric, or the axial vector might 
have more than two poles. In such extreme cases, this 
would not work. But we shall see through our numerical 
simulations that there is a significant dynamical regime 



where exactly 2 poles exist, and these problems do not 
arise. 

Given the location of the two poles on the horizon, how 
do we assign a 3-vector to them, and thereby obtain all 
the components of the spin vector S? An obvious start- 
ing point would be to use the unit normal vector R a at 
the poles. This would not give a unique answer in the 
absence of reflection symmetry. Alternatively, we could 
consider the curl of ip a suitably averaged over the hori- 
zon. However, even if we could successfully assign such a 
3-vector uniquely, it is not clear in general how this vec- 
tor should be compared with the spin 3-vector calculated 
at spatial infinity. This could be done in spacetimes with 
global axisymmetry, but this is not available to us in the 
present case. In the absence of a solution to this problem, 
we simply define the direction of the spin to be the Eu- 
clidean unit-norm vector tangent to the coordinate line 
joining the two poles. The spin- vector Sm is then equal 
to this Euclidean unit-norm vector multiplied by the Iso- 
lated Horizon spin obtained from Eq. The definition 
of S might need to be further refined, however it seems 
to be satisfactory for our purposes. This definition of the 
spin vector reproduces the Bowen-York spin parameters 
on the initial slice, and should remain reasonable as long 
as the coordinates do not become too distorted. In addi- 
tion to using the Killing vector ip a , we also found it useful 
to define angular momenta with the flat space coordinate 
rotational killing vectors 

= [°>-(z- z c),(y ~ Vc)] , 
¥>% = [{*- z c ),0,-(x - x c )] , 
V% = {-{y-Vc),{x-x c ),Q), 

where (x c ,y c ,z c ) is the coordinate centroid of the hori- 
zon. We can then obtain the coordinate-base spin vec- 
tor ^coord = (S x , S y , S z ) by replacing the approximate 
Killing vector in Eq. (J4| with the three coordinate ro- 
tational vectors (i.e. Si = ± § AH (tpfR b K ab ) d 2 V). This 
definition of the spin direction reproduces the Bowen- 
York spin parameters on the initial slice as well, and pro- 
duces reasonable results at later times for the gauges used 
here. (Of course this latter coordinate based calculation 
will not yield an accurate evaluation of the spin direction 
or magnitude for more general gauges, while the former 
approximate Killing vector calculation will produce accu- 
rate spin magnitudes for generic gauges.) In both cases 
the spin magnitude is the Euclidean norm of the 3-vector 
S. Note that in the former case this Euclidean norm is 
precisely the spin given by Eq. (|4|). 

We solve for the approximate Killing vector field ip a on 
the horizon using standard spherical-polar coordinates. 
In these coordinates the Killing vector is obtained with 
highest accuracy when its poles are aligned with the co- 
ordinate poles of the (9, 0) coordinates. To make the 
calculation as accurate as possible, we find the minimum 
of (p a (f a in the northern hemisphere and then rotate the 
angular coordinates so that the new north pole is aligned 
with the minimum of <p a ip a - We then re-calculate (p a to 
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TABLE I: Initial data for quasi-circular, equal-mass black-hole binaries. The binaries have an ADM mass of (1.0000±0.0005)M, 
with orbital frequency Mfl fixed to 0.0500, and initial proper separations I. The punctures are located at (±Jf, 0, 0), with mass 
parameter m p , momentum (0, ±P, 0), spin angular momentum (0, S y , S z ), and specific spin S/m 2 (m is the horizon mass). 



Name 


Sy/M 2 


S z /M 2 


X/M 


P/M 


J/M 2 


S/m 2 


l/M 


m p /M 


SP3 


0.128725 





3.276347 


0.133587 


0.91243 


0.5013 


10.20 


0.43025 


SP4 


0.091198 


0.091198 


3.179908 


0.1314406 


1.03454 


0.5007 


9.94 


0.43037 


SO 


0.0 


0.0 


3.280 


0.1336 


0.876 


0.0 


10.01 


0.4848 


SC 


0.0 


0.025757 


3.2534 


0.1330 


0.917 


0.1001 


9.93 


0.4831 


S++ 


0.0 


0.1939 


3.0595 


0.1291 


1.1778 


0.757 


9.27 


0.3344 


s- 


0.0 


-0.1924 


3.465 


0.1382 


0.5729 


-0.757 


10.3 


0.3344 



obtain a more accurate location of the minimum and it- 
erate until the new minimum of ip a ip a lies on the north 
pole. There is a complication in the above procedure 
in that we cannot calculate ip a on the coordinate poles 
themselves (since the 2D Christoffel symbol is singular). 
In practice we stop iterating when the minimum of ip a ip a 
lies within 2 angular grid-points of the coordinate pole. 
The spin-direction associated with the minimum of tp a ip a 
therefore cannot be obtained with higher precision than 
a few angular grid sizes. It might be possible to im- 
prove the accuracy by considering multiple patches on 
the horizon to avoid the coordinate singularity, or to use 
a spectral decomposition. 

We found that using 160 points in the 9 direction and 
320 points in the <fi direction provides reasonable results 
for the horizon spin calculation, with errors in the spin 
direction of about 2°. Adding significantly more points 
only increases the numerical error because the horizon 
algorithm itself uses far fewer points to locate the appar- 
ent horizons and the underlying numerical grid has a far 
coarser resolution. 

The configurations discussed in this paper contain ei- 
ther Pl-symmetry, i.e. (x,y,z) — > (— x,— y, z), or parity- 
symmetry, i.e. (x, y, z) — > (—x, —y, —z). We exploit these 
two symmetries to reduce the grid size by a factor of 2. 
The zero-spin and (anti-)aligned spin binaries have the 
additional symmetry (x,y,z) — ► (x, y, — z). We imple- 
ment the parity-symmetry boundary conditions using a 
locally modified version of the Pi-symmetry boundary 
thorn kindly provided to us by Erik Schnetter. Note 
that behavior of the components of g a b and k a {> un- 
der these symmetries can be obtained from the behav- 
ior of these components under simple reflections. Thus, 
g xz (x 7 y,z) = g xz (-x 1 -y,-z) under parity-symmetry 
while g xz (x 1 y, z) = -g xz (-x, -y, z) under Pi-symmetry. 
The Einstein equations preserve these symmetries. In the 
case of spinning, equal-mass binaries, parity-symmetry 
requires that the two spins be equal in magnitude and 



parallel (i.e. the spin- vector behaves as a pseudo vector). 



IV. INITIAL CONFIGURATIONS 

We study two configurations of non-aligned-spin bina- 
ries with parallel spins (equal in magnitude) that exhibit 
spin and orbital-plane precession, as well as spin-flips of 
the remnant horizon spin with respect to the individual 
horizon spins. We choose configurations where the binary 
separation is small enough that the spin-orbit coupling is 
large, but large enough that the binaries complete at least 
~ 1.75 orbits before merging. The first configuration, 
which we denote with SP3 starts with the spins aligned 
along the initial orbital plane. This can be interpreted as 
a binary in which one black hole orbits about the pole of 
the second black hole. The second configuration, which 
we denote with SP4, starts with the spins pointing 45° 
above the initial orbital plane, corresponding to a binary 
in which infall occurs initially along a plane tilted with 
respect to both spins. In both cases the masses and spins 
of the two holes are equal (i.e. spins parallel and equal 
in magnitude). Setting the two masses and spins equal 
ensures that the system is parity-symmetric, but still is 
generic enough to display both spin and orbital plane 
precession as well as a spin-flip in the direction of the 
orbital angular momentum. The initial data parameters 
for these two configurations, which were obtained using 
the 3PN equations of motion, are given in Table Q] We 
also report the initial-data parameters for the previously 
studied aligned (S++, SC), anti-aligned (S--), and non- 
spinning binaries (SO). The PN data provides the punc- 
ture location, momenta, and spins. We complete the data 
by choosing puncture mass parameters (equal for the two 
punctures) such that the total ADM mass of the system 
is 1. 



V. RESULTS 



We evolved the SP3 configuration using central res- 
olutions of h = M/22.5, h = M/25, and h = M/30; 



with grid-sizes of 576 2 x 288, 640 2 x 320, and 768 2 x 384 
respectively. We used 'multiple transition' fisheye trans- 
formation [30( to place the outer boundaries at 200Af; 
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far enough away that boundary effects do not interfere 
with the orbital dynamics of the system. In addition, we 
also evolved the SP4 configuration with a central resolu- 
tion of h = M/25, a grid size of 640 2 x 320, and outer 
boundary at 200 Af. 

Figures Q] and [5] show the puncture trajectory and 
horizon-spin direction along this track for the SP3 con- 
figuration (the latter suppressing the z-direction). Note 
that the scale of the z-axis in Fig. [T] is l/10th that of 
the x and y axes. From the plots one can clearly see the 
orbital plane precess out of the equatorial plane, as well 
as the spin axis rotating by approximately 90° in the xy 
plane during the course of the merger. The spins are ini- 
tially aligned along the y-axis, but at merger they show 
both a significant z-component and an approximate 90° 
rotation to the — x-axis. The individual horizon spins 
at the merger are S coord = (-0.121 ± 0.002, -0.007 ± 
0.003, 0.037 ± 0.003) (we use the coordinate based mea- 
sure of the spin at the merger because the calculation of 
Sih is not accurate when the black holes are this close 
together; see comments below). Hence the total preces- 
sion angle for the SP3 configuration is 8 P = 98°. Note 
that there is no discernible correlation between the ori- 
entation of the projected horizon and the projected spin 
direction. 

In Fig. [3] we show the coordinate (S COO rd) and Killing 
vector based (Sih) calculation of the spin components 
versus time. Sih displays a step-function-like behavior 
due to the difficulty in finding the poles (i.e. the ze- 
roes of (fi a fa) in the Killing vector accurately. As dis- 
cussed above, the Killing vector calculation is most ac- 
curate when its poles are located at the coordinate poles 
of the (6, </>) coordinates on the horizon. However, the 
difficulty in calculating the Killing vector itself near the 
coordinate poles introduces an uncertainty in the loca- 
tion of the Killing vector poles. We also stress that the 
direction associated with the location of the two poles 
of the Killing vector is coordinate dependent. From the 
figure we see that the cc-and-y-components of Sih oscil- 
lates about the much more regular x and y components 
of Dcoord) while the z-component of 5m is consistently 
larger than the z-component of S COO rd- The calculation 
of <Sm (and hence Sih) breaks down prior to the merger 
when the horizons get too close (and hence the mutual 
tidal distortions destroy the approximate axial symme- 
try). Scoord, however, continues to produce reasonable 
results through the merger. Thus, it is S coor d that shows 
the clear rotation of the spin from the y-axis to the x- 
axis at the merger. Note that the uncertainties in the 
spin directions do not correspond to uncertainties in the 
magnitudes of the spin. For the Killing vector based 
calculation Sih, it is the spin magnitude that is deter- 
mined with high accuracy. Figure |4] shows between third 
and fourth-order convergence of the components of S coor d 
from the three resolutions (the third-order error may be 
due to third-order errors leaking out of the puncture as 
well as third-order errors from the horizon calculations), 




FIG. 1: The puncture trajectories along with spin direction 
(every 4M) for the SP3 configuration for the Af/30 resolution 
run. The spins are initially aligned along the y-axis, but rotate 
by ~ 90° during the 1.25 last orbits and also acquire a non- 
negligible z-component. Note that the z-scale is l/10th the x 
and y scale. 



while Fig. [5] shows the value of the z-component of the 
specific spin S z /m 2 (where m is the horizon mass) based 
on the z-component of Sih for the three resolutions. In 
this latter figure the curves have been translated. (A 
convergence plot of Sih would not be meaningful be- 
cause the size of the step discontinuities in Sih are larger 
than the differences in the spin direction with resolution.) 
For the z-component of the spin, we expect that, given 
the lack of significant oscillations plaguing the x and y 
components, Sm gives a better measurement than the 
more highly coordinate-dependent S COO rd- The spin-up 
of the z-component of the specific spin by 0.16 is ten 
times larger than the analogous spin-ups of about 0.01 
seen in the zero-spin (SO) and aligned-spin (SC) config- 
urations (see [34|]). However, as can be seen in Fig. [3l 
the spin magnitude does not increase significantly. Thus 
this spin-up in the z-direction is not equivalent to the 
spin-up observed in the case of the aligned-spin and non- 
spinning binaries. In those cases the spin-up involved an 
increase in the spin magnitude, while here it primarily 
involves a rotation of the spin vector out of the xy plane. 
This rotation of the spin out of the xy plane follows the 
post-Newtonian predictions of Eq. ([TJ) . 

In Fig. [6] we plot the x and y components of the spin 
as a function of the z component for both the Post- 
Newtonian predicted spins (using numerical tracks) and 
the numerically determined spins. Plotting the data in 
this manner removes the ambiguity of assigning the ap- 
propriate Post-Newtonian time to the numerical time co- 
ordinate on the horizon. The qualitative behavior of the 
spin in our numerical simulation is consistent with the 
post Newtonian spin for most of the evolution (smaller 
values of S 2 ). 

The puncture trajectories are third-order convergent as 
is demonstrated in Fig. [7J The x-component of the track 
appears to show poorer convergence between t — lOAf 
and t — 6QM but this is likely due to the coarseness 
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FIG. 2: The projection of puncture trajectories and spin for 
the SP3 configuration onto the xy plane along with the indi- 
vidual apparent horizons for the M/30 run. The horizons and 
spins are shown at t = 0, 20M, 40M, 160M, 164M. The 
first common horizon (also shown) formed at t — 164. 2M. 
The spins are initially aligned along the j/-axis but rotate by 
~ 90° during the last 1.25 orbits. The spin of the second 
black hole (not shown) is equal to the first. 
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FIG. 3: The spin components and magnitude versus time 
for the SP3 configuration as calculated using the coordinate 
rotational vectors (coord) and the poles of the approximate 
Killing vector (IH) for the M/30 resolution. The calcula- 
tion of the approximate Killing vector breaks down near the 
merger (which occurs at t — 164.2), but the purely coordi- 
nate based calculation continues to produce reasonable re- 
sults. Note that the direction obtained from the Killing vec- 
tor oscillates about the coordinate based direction. Also note 
the spin has just rotated by 90° in the xy plane at the time 
of merger. The spin magnitude remains essentially constant 
throughout the merger phase. The magnitude of the spin cal- 
culated from the Killing vector is coordinate invariant and, 
unlike the spin direction, is expected to be more accurate 
than the coordinate based calculation. 



of the grid. Note that the trajectories are calculated 
by integrating dtx a = —j3 a [22] at the location of the 
punctures. The curve is expected to converge to lower 
order because the shift is not smooth on the puncture 
(see [lij for a discussion on the behavior of the evolved 
fields at the punctures). However, as is shown below, the 
constraint violations also show third-order convergence, 
which indicates that lower-order convergent effects leak 
out of the puncture. These lower-order errors are likely 
not observed in the waveform (see below) because of a 
larger fourth-order error term dominating the third-order 
error terms at these resolutions. 

The final remnant horizon for the SP3 configuration 
has mass M n /M = 0.9613 ± 0.0007 with specific spin 
S/M^ = 0.7215 ± 0.0003. The remnant spin compo- 
nents calculated from the approximate Killing vector are 
Sm/M 2 = (-0.045 ± 0.001, 0.199 ± 0.003, 0.638 ± 0.003) 
with magnitude Sm/M 2 = 0.669±0.001. The total ADM 
mass and angular momentum (i.e. initial mass and angu- 
lar momentum) of the system are Madm/M = 1.00000± 
0.00005 and J AI)M /M 2 = (0,0.257450,0.875352). Hence 
(3.87 ± 0.07)% of the mass and (23.6 ± 0.1)% of the 
angular momentum were converted into radiation, with 
SJ Z /M 2 = (0.237± 0.003). The system gained net angu- 
lar momentum in the x-direction but lost (22.7 ± 0.4)% 
of its angular momentum in the y-direction and (27.1 ± 
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FIG. 4: A convergence plot of the coordinate-base S CO ord cal- 
culation for the SP3 configuration. Note that for our choice 
of resolutions, third-order convergence is demonstrated by 
5(M/22.5) - S{M/2b) = (S(M/25) - S(M/30))C, where 
C = 0.88. The spin is initially third-order convergent, with 
higher order-convergence apparent at later times. 



8 



0.18 
0.16 
0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 




M/22.5 
M/25 
M/30 




160 



FIG. 5: The z-component of the Killing vector based calcula- 
tion of the spin rescaled by the square of the horizon mass for 
the three resolutions. The curves have been translated by a 
distance equal to the difference in merger times of the Af/22.5 
and M/25 runs with the merger time of the Af/30 run. In this 
configuration, unlike the previously studied aligned-spin bi- 
nary, the 'spin-up' appears to be large. However, in this case 
the 'spin-up' is actually the rotation of the nearly-constant 
spin-vector towards the z-axis, rather than an increase in the 
spin magnitude. 
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ponents of the spin of the individual horizons for the SP3 
configuration as functions of the S z component of the spin 
(which increases monotonically in time prior to the merger). 
Note the very reasonable agreement for most of the evolution. 
The plot terminates at the merger. 



0.3)% of its angular momentum in the z-direction. Thus 
the binary preferentially radiated angular momentum 
in the direction of the initial orbital angular momen- 
tum. We also obtained estimates of the radiated mass 
and z-component of the angular momentum from ip4 of 
(3.8 ± 0.1)% and (0.24 ± 0.02)M 2 respectively, in exce- 
lent agreement with the results from the remnant horizon 
parameters. (The errors in the radiated mass and angu- 
lar momentum from the waveform are relatively large 
due to boundary reflection contaminating the late-time 
waveform.) Note that the excellent agreement in SJ Z be- 
tween the final horizon direction measurement and the 
radiated z-component of the angular momentum indi- 
cates that we obtain the final horizon spin magnitude 
and direction to within the expect 2° accuracy. For com- 
parison we also give the remnant spin direction calcu- 
lated using the coordinate rotation vectors S coor d/M 2 = 
(-0.033±0.002, 0.190±0.001, 0.6395±0.0003) (where the 
errors are a measure of the flatness of the components of 
Scoord versus time). The corresponding spin magnitude 
S COOI d/M 2 = 0.668 ±0.002 agrees with the norm Sm and 
the directions agree to within 1.3° (the expected error in 
the direction determination is ~ 2°). 

We measure the angle of the spin-flip both with respect 
to the initial individual horizon spins and the individual 
horizon spins at the merger. For the SP3 configurations 
these angles are 72° and 71° respectively. In Fig. [5] we 
show the spin direction of the individual horizons and 
spin direction of the remnant horizon. The smooth pre- 
cession and discontinuous flip are apparent. Note that 
the spin flip, unlike the spin precession, cannot be mod- 
eled accurately by a post Newtonian expansion due to the 
highly non-linear merger process that converts roughly 
25% of the initial total angular momentum into gravi- 
tational radiation. The final spin direction for the SP3 
configuration is rotated by 4.0° with respect to the initial 
ADM angular momentum. This rotation, though small, 
is larger than the expected error in our spin-direction al- 
gorithm, indicating that there is a small net change in 
the direction of the angular momentum. 

The spin direction and puncture trajectories are coor- 
dinate dependent measures of precession. The waveform, 
on the other hand, should provide a coordinate indepen- 
dent measure of the precession. To show the effect of pre- 
cession on the waveform we examine the (£ = 2, to = 1) 
mode. This mode vanishes identically for the zero-spin, 
aligned-spin, and anti-aligned spin cases previously stud- 
ied (i.e. SO, SC, S++, S--) 0, [H H HI- In Fig. M 
we show the [l = 2, m = 1) mode at an extraction ra- 
dius of r = 15M as well as a convergence plot of this 
mode (showing fourth-order convergence). The contri- 
butions of the (£ — 2,m = ±1) modes to the radiated 
mass are smaller than the contributions of the dominant 
(I = 2, to = ±2) modes by a factor of ~ 20, while the con- 
tributions of the (£ = 2, to = ±1) modes to the radiated 
angular momentum are smaller than the contributions of 
the dominant (£ = 2, m = ±2) modes by a factor of ~ 80. 

We next examine how the results change when we set 
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FIG. 7: The difference in track locations between the M/22.5 
and M/25 runs as well as the difference in track locations 
between the M/25 and M/30 runs. The latter difference has 
been rescaled by C = 0.88 to demonstrate third-order con- 
vergence. The reduced order of convergence for X between 
t = 10M and t = 60M is likely due to the coarseness of the 
grid. 
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FIG. 9: The real and imaginary parts of the (I = 2, m = 1) 
mode of ip4 for the SP3 configuration with resolution h = 
M/25 (upper panel), as well as a convergence plot of this 
waveform (lower panel). Note that the waveform is fourth- 
order convergent (as is evident by the agreement of the two 
differences) . 




FIG. 8: The spin-direction of the individual horizons every 
4M during the spin-precession phase and the final horizon 
spin-direction for the SP3 configuration. The arrows indi- 
cate the spin-direction only, not the magnitude. Note the 
continuous change in the spin-direction during the precession 
stage and the discontinuous jump (or flip) to the remnant 
spin-direction. 



the initial spins closer to the z-axis. The SP4 configu- 
ration has the same total spin as the SP3 configuration 
but at an angle of 45° with respect to the orbital plane. 
This rotation of the spin has two significant effects. First, 
the SP4 configuration has a significant spin in the same 
direction as the orbital ang ular momentum, and from 
our previous results [34], |35| | we expect that the binary- 
merger will be delayed due to the resulting spin-orbit re- 
pulsive effect. Second, the amplitude of the orbital plane 
precession will be reduced (i.e. there is no orbital plane 
precession if the spins are rotated 90° with respect to 
the orbital plane, and the amount of precession should 
vary smoothly with angle). In Fig. [TU] we show the z- 
component of the first puncture trajectory versus time 
for SP3 and SP4, where the latter has been rescaled by 
V2. Note that at early times the rescaled tracks agree 
perfectly. Thus, for a given magnitude of the spin, the 
orbital plane precession has a sin$ dependence, where 
•d is the angle between the spin and orbital angular mo- 
mentum. At later times, the spin-orbit coupling induced 
delay in the merger becomes evident and the two tracks 
no longer agree. 

Comparing the waveforms from the SP3 and SP4 con- 
figurations is complicated by the fact that the initial 
data burst masks the early-time behavior. Neverthe- 
less, there is a small region just after the initial pulse 
leaves the system (see Fig. ITTj) where it is evident that 
the (£ = 2, m = 1) mode scales with sin$. However, 
at later times the differences in the orbital dynamics de- 
stroys this scaling. 

In Figures Q2J [131 [14] we show the spin magnitude and 
direction for the SP4 configuration as well the 3D punc- 
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FIG. 10: The ^-component of the puncture trajectory (for 
the puncture initially located at x > 0) versus time for the 
SP3 and SP4 configurations with resolution h — M/25. Note 
that after rescaling by y/2 the two trajectories agree at early 
times. At later times the spin-orbit coupling in the partially 
aligned case delays the merger; causing the two trajectories 
to diverge. Both curves terminate at their respective merger 
times. 



FIG. 12: The spin components and magnitude of the indi- 
vidual horizons (the two horizons have equal spins) for the 
SP4 configuration up to merger (with resolution h — M/25). 
Note that in this configuration the spins rotate by 135° in 
the xy plane prior to merger. Also note that the significant 
z-direction spin-up is once again caused by a rotation further 
out of the xy plane. Once again the Killing vector determi- 
nation for the spin direction oscillates about the coordinate 
definition for the spin direction. The approximate Killing 
vector calculation begins breaks down at around t — 180 M . 
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FIG. 11: The imaginary part of the (£ = 2, m = 1) compo- 
nent of ip4 for the SP3 and SP4 configurations extracted at 
r — 15M with a central resolution of M/25. Note that be- 
tween t = 50 and t — 75 the two waveforms scale with sin # as 
is evident by the good agreement between the two waveform 
after rescaling the SP4 waveform by v2 (i.e. 1/ sin(7r/4)). 
However, this scaling breaks down at later times (t > 85Af) 
due to the differences in orbital decay arising from the in- 
creased stability of the SP4 configuration. This scaling also 
breaks down at early times because of the non-physical initial 
data radiation pulse. 



ture trajectories and spins and a projection of the trajec- 
tories and spin direction onto the xy plane. Due to the 
increased stability of aligned spin binaries, this config- 
uration completes 2| orbits prior to merger (compared 
to if for SP3). Consequently spin-precession rotates the 
spin vector by an additional 45° compared to the SP3 
configuration. In this case the Killing vector based cal- 
culation of the spin remains accurate long enough that 
the spin-precession rotation beyond 90° (i.e. the local 
minimum in S x is observed in Sm as well as iScoord)- 
Note that, once again, there is a significant spin- up in 
the z-direction caused by a rotation of the spin vector 
towards the z-axis (rather than a net increase in the spin 
amplitude), and that there is no discernible correlation 
between the projected horizon orientation (i.e. the ori- 
entation of the semi-major axis) and the projected spin 
direction. The individual horizon spins at the merger are 
Scoord = (-0.033 ± 0.005, -0.041 ± 0.003, 0.114± 0.001). 
Hence the total precession angle for the SP3 configura- 
tion is 6 P = 151°. 

The remnant horizon for the SP4 run formed just as the 
boundary reflections began to contaminate the interior. 
Consequently the error bounds for the mass and spin of 
the remnant are higher for SP4 than SP3. The final rem- 
nant mass is M H = (0.9524 ± 0.0002)M with a spin pa- 
rameter of S/Af£ = 0.805 ± 0.002. The spin components 
are S comd /M 2 = (-0.020 ± 0.003, 0.121 ± 0.002, 0.720 ± 
0.002), where we used the coordinate-base definition to 
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FIG. 13: The puncture trajectories and horizon spin direction 
(shown every AM until merger) for the SP4 configuration. 
Note that the 2-scale is l/10th that of the x and y scales. 
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FIG. 14: The xy plane projection of the orbital track, 
apparent horizons, and spin direction for the SP4 config- 
uration up to merger. The horizons are given at t = 
0, 20M, AOM, 180M, 196M. The common horizon formed 
at t — 195. AM. The projected spin vector decreases in mag- 
nitude at late time due to the spin rotating further out of the 
xy plane. Note that the spin direction precesses by 135° in 
the xy plane during the merger. The spin of the second black 
hole (not shown) is equal to the first. 



calculate S. (The Killing vector based calculation of the 
spin could not be obtained accurately because the system 
became approximately axisymmetric after the boundary 
errors affected the remnant spin parameters.) The ini- 
tial ADM mass and angular momentum for this system 
were M ADM /M = 1.00000 ± 0.00005 and J AD m/M 2 = 
(0,0.182396,1.01833). Hence, (4.76 ±0.02)% of the mass 
and (29.4 ± 0.2)% of the angular momentum were con- 




FIG. 15: The spin-direction of the individual horizons every 
AM during the spin-precession phase and the final horizon 
spin-direction for the SP4 configuration. The arrows indi- 
cate the spin-direction only, not the magnitude. Note the 
continuous change in the spin-direction during the precession 
stage and the discontinuous jump (or flip) to the remnant 
spin-direction. 



verted into gravitational radiation. The system gained 
net angular momentum in the —x direction, while los- 
ing (34 ± 1)% and (29.3 ± 0.2)% of its angular momen- 
tum in the y and z directions respectively. It thus ap- 
pears that this configuration preferentially radiates angu- 
lar momentum in the orbital plane. However, we caution 
the reader that the errors quoted for the SP4 configu- 
ration for the mass and spin do not take into account 
either possible boundary effects or finite-difference trun- 
cation errors. The effect of the radiation on the angular 
momentum direction is small, with the SP4 spin direc- 
tion rotated by only 1.5° with respect to the initial ADM 
angular momentum. 

The spin-flip angle of the remnant spin with respect 
to the initial individual spins and the individual spins at 
merger are 35° and 32° respectively. In Fig. [15] we show 
the spin-direction of the individual horizons and spin- 
direction of the remnant horizon. The smooth precession 
and discontinuous flip are apparent. Note that in this 
case boundary reflections contaminate the waveform at 
large r prior to the merger. Consequently we do not 
obtain reliable measurements for the radiated mass and 
angular momentum from the waveform. 

As was mentioned above, the calculation of the spin 
direction is coordinate dependent. Nevertheless, these 
particular coordinates show remarkable agreement be- 
tween the puncture trajectories and the waveform. In 
Fig. [16] we show the orbital part of the (t = 2,m = 2) 
component of ip4 extracted at r = 10 M , where we trans- 
lated the SP4 waveform and multiplied by a constant 
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phase in such a way that plunge part of the waveforms 
agree [32]. Note that the SP3 configuration shows ~ 4 
cycles of orbital motion prior to the plunge (i.e. the last 
trough in the plot), while the SP4 configuration shows 
~ 5 cycles. Thus we expect that the SP4 track should 
contain approximately one-half of an orbit more than the 
SP3 configuration. Although there is an uncertainty in 
the exact location of the beginning of the 'plunge' wave- 
form in the figure, its approximate location will be given 
by the formation time of the common horizon at this 
resolution plus the coordinate distance to the extraction 
sphere (here we identify the start of the plunge with the 
first trough located at t = 168 = T C ah + 10M), the 
number of cycles in both configurations prior to the last 
peak shown is consistent with 1 /4 of an orbit more than 
the number of orbits (i.e. 1.75 and 2.25 respectively) ob- 
served in Figs. [2] and 1141 Interestingly, the number of 
cycles after the initial pulse of radiation (3.5 and 4.5 re- 
spectively) is in excellent agreement with the number of 
orbits observed in the puncture trajectories. Thus these 
coordinates appear to reasonably reproduce the orbital 
dynamics of the binary. This fidelity by which the co- 
ordinates reproduce the merger dynamics, and the rel- 
atively good agreement between the two measurements 
of the spin direction, is the motivation for using these 
coordinate dependent measurements to measure the spin 
direction. In addition, the very good agreement for the 
radiated z-component of the angular momentum based 
on the difference between the remnant Sm and Jadm for 
SP3 and the waveform-based calculation, indicates that 
this method provides an accurate measurement of the 
remnant spin direction and magnitude. We expect the 
method to provide more accurate results for the remnant 
horizon than the individual horizon because the remnant 
spacetime is axisymmetric, and, as pointed out above, 
there is a natural way to assign a spin direction to hori- 
zons when the spacetime is axisymmetric. 

There appears to be a small trend towards spin- 
angular momentum alignment. To observe this effect we 
measure the angle between the spin and final angular mo- 
mentum (i.e. the remnant spin), 9js = arccos(5 l • Jf). In 
Fig. [17] we show 6js versus time for the SP3 and SP4 
configurations using both Sm and Scoord- Both of these 
measurements show increasingly larger oscillations in Ojs 
versus time. The Sm based measurements indicate that 
the spins in the SP3 configuration migrate ~ 6° towards 
Jf , while the spins in the SP4 configuration migrate ~ 8° 
towards Jf (we measure the angle at the approximate 
midpoint of the oscillation). However, the iScoord based 
measurements indicate that the spins in the SP3 con- 
figuration migrate by only ~ 3°, while the spins in the 
SP4 configuration migrate by only ~ 2.5°. It is unclear 
which measurement is superior. On the one hand S^oord 
is more strongly coordinate dependent, but, on the other 
hand, the horizons become more distorted as they ap- 
proach each other, which increases the error in Sjjj (i.e. 
the horizons deviate increasingly strongly from axisym- 
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FIG. 16: The real part of the pre-merger (I = 2, m = 2) 
component of ipi extracted at r = 10M. The SP4 waveform 
has been translated by 37M and multiplied by a constant 
phase factor prior to taking the real part. The plunge part 
of the waveform begins roughly at t = 168M. Thus there are 
3.5 cycles of orbital radiation prior to the plunge (but after 
the initial data pulse) for SP3 and 4.5 for SP4. This number 
of cycles match the number of orbits in Figs. [2land [T4l The 
initial data pulse appears to mask an additional 1/2 cycle of 
orbital motion. 



TABLE II: The angle between the initial spin direction and 
orbital angular momentum §, the total spin precession an- 
gle O p , spin flip angle between the initial spin direction and 
remnant spin direction Qfa P , and the net change in the z- 
component of the angular momentum as calculated using the 
Isolated Horizon spin direction of the remnant <5J Z (IH) and 
waveform 5Jz{^i) for the SP3 and SP4 configurations. 



Config 


■d 


Q P 


Oflip 


sj x (m) 


5Jz(i>4.) 


SP3 


90° 


98° ± 2° 


72° ± 2° 


0.237 ± 0.003 


0.24 ±0.02 


SP4 


45° 


151° ±2° 


35° ± 2° 







metry). This is an interesting effect which we plan to 
study in much more detail with improved techniques to 
measure the spin direction. We summarize the main re- 
sults from spin-direction calculations for the SP3 and SP4 
simulation in Table HT1 

The stability of the spinning binaries is strongly depen- 
dent on the direction of the spin. In Table IIIII we show 
the merger times (Tqah) of the SP3 configuration versus 
resolution and an extrapolation to infinite resolution, as 
well as the merger time for the single SP4 run. The 'ex- 
trapolated' value of the SP4 merger time was computed 
by adding the difference between the extrapolated and 
h = A//25 merger times for the SP3 configuration to the 
h = M/25 merger time for the SP4 configuration. The 
extrapolated values of 176±3 and ~ 214, for SP3 and SP4 
respectively, are in large part consistent with the results 
from the aligned-spin binaries if we replace S in Eq. (|A1[) 
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FIG. 17: The angle in degrees between the individual appar- 
ent horizon spins and the final total angular momentum for 
the SP3 and SP4 configuration. The SP3 angles have been 
translated by —37.21° so that the translated SP3 and SP4 an- 
gles agree at t = 0. The curves terminate at the point when 
the Isolated Horizon algorithm becomes inaccurate. 
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FIG. 18: The convergence of the Hamiltonian constraint vio- 
lation for the SP3 configuration along the a>axis at t = 76M 
when the punctures cross the x-axis for the second time. 
The Hamiltonian constraint shows third-order convergence. 
Points inside the domain of dependence of the boundary have 
been excluded. The high-frequency features near the outer 
boundary are due to the extreme fisheye de-resolution and 
converge with resolution. 



TABLE III: Merger times for the SP3 and SP4 configuration 
versus resolution, as well as an extrapolation to infinite res- 
olution (see text for an explanation of the 'extrapolation' of 
the SP4 result.) 



Resolution 


SP3 


SP4 


Af/22.5 


152.0 ±0.2 




M/25 


157.4 ± 0.2 


195.4 ±0.2 


M/30 


164.2 ±0.2 




M/oo 


176 ±3 


~ 214 



with 2S z /m 2 (S is the total spin, hence the factor of 2). 
The predicted merger times are Tcah = 172 ± 1 and 
T C ah = 200 ± 2 for the SP3 and SP4 configurations re- 
spectively. The differences between these predictions and 
the actual extrapolated merger times can be explained by 
the net rotation of the component spins towards the z- 
axis, which helps stabilize the binaries. 

We conclude this section by showing that the con- 
straint violations converge to third-order. Although the 
code uses purely fourth-order stencils, lower order errors 
both from the lower differentiability of the evolved fields 
at the punctures, as well as from the second-order ac- 
curacy of the initial data, lead to a global third-order 
error in the constraint violation. Figures fl~8l \19\ l20l show 
the Hamiltonian constraint, Momentum constraint, and 
BSSN constraint (G l = f l +dfg li ) violations at t = 76M 
(the time when the punctures cross the x-axis for the 
second time) along the x-axis for the SP3 configuration. 
The constraint violations have been multiplied by (hi/h) 3 
(where hi — Af/22.5) to demonstrate third-order con- 



vergence. Note that the high-frequency features near 
the outer boundary are due to the extreme fisheye de- 
resolution near the outer boundary and converges to zero 
with resolution. 



VI. CONCLUSION 

In this paper, we evolved systems of equal-mass and 
equal-spin black-hole binaries with initial spins aligned 
perpendicular to, and 45° to, the orbital angular mo- 
mentum. We observed the combined effects of spin and 
orbital plane precession as predicted by post Newtonian 
theory, with dramatic, total precessions of ~ 98° and 
~ 151° in the SP3 and SP4 simulations respectively. 
Both configurations resulted in large spin flips between 
the individual horizon spin directions and the final rem- 
nant direction, with the SP3 configuration showing a spin 
flip of ~ 72° and the SP4 configuration (which had spins 
initially more closely aligned with the orbital angular mo- 
mentum) resulting in a spin flip of ~ 34°. We see a 
possible small trend to spin — orbital angular momentum 
alignment, but no evidence for significant spin- up of the 
black-hole spins. 

Although the configurations studied here are parity- 
symmetric, this symmetry was only chosen to reduce the 
memory footprint of the simulations (allowing for higher 
resolution runs); it does not affect the stability of the 
'moving punctures' algorithm. Notably, the punctures 
move out of the xy plane and can get arbitrarily close 
to the numerical gridpoints. Despite this symmetry, the 
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FIG. 19: The convergence of the Momentum constraint vio- 
lation for the SP3 configuration along the a;-axis at t = 76M 
when the punctures cross the x-axis for the second time. The 
Momentum constraint shows third-order convergence. Points 
inside the domain of dependence of the boundary have been 
excluded. In each panel the solid (black) curve, dotted (red) 
curve and dashed (blue) curves are the M/22.5, M/25 and 
A//30 constraint violations respectively. The constraints have 
been rescaled by (hi/h) 3 (hi — M/22.5). The high-frequency 
features near the outer boundary are due to the extreme fish- 
eye de-resolution and converge with resolution. 



FIG. 20: The convergence of the BSSN constraint violation 
for the SP3 configuration along the a;-axis at t — 76M when 
the punctures cross the :r-axis for the second time. The BSSN 
constraint shows third-order convergence. Points inside the 
domain of dependence of the boundary have been excluded. 
In each panel the solid (black) curve, dotted (red) curve and 
dashed (blue) curves are the M/22.5, M/25 and M/30 con- 
straint violations respectively. The constraints have been 
rescaled by (hi/h) 3 (hi = M/22.5). The high-frequency fea- 
tures near the outer boundary are due to the extreme fisheye 
de-resolution and converge with resolution. 



SP3 and SP4 configurations display most of the signifi- 
cant spin-orbit coupling effects associated with spinning 
binaries: spin and orbital plane precession, spin flips, and 
enhanced stability of the semi-aligned configuration. The 
only significant spin-orbit coupling effect not shown by 
these configurations is a spin-orbit induced kick of the 
remnant hole. In order to see these kicks we would need 
to evolve configurations without parity symmetry. 

Our methods for calculating the spin direction produce 
reasonable results for our choice of gauge conditions. Fu- 
ture work will concentrate on improving this calculation 
with alternative choices of the gauge parameters (e.g. r\ 
in the Gamma-driver shift and various different choices 
of initial values for the lapse function) and with alterna- 
tive forms of the Gamma-driver shift condition. Notably, 
the lack of agreement between the location of the hori- 
zon semi-minor axis and the spin direction indicates that 
these coordinates are not yet ideal. In addition, it would 
be useful to calculate independent measures of the qual- 
ity of the approximate Killing vector, for example the 
norm of the Lie derivative of the 2-metric on the horizon 
£tpqab- It is interesting to note that the purely coordi- 
nate measurement S'coord gives reasonable results for the 
spin direction and amplitude, and since this calculation 
is both more robust (i.e. the approximate Killing vector 
may not exist) than the approximate Killing vector cal- 



culations and easier to implement, it may prove to be 
a convenient measurement of the spin for those codes 
that have not implemented the approximate Killing vec- 
tor finding algorithm. 

From the mathematical side, further investigation is 
required to make the definition of the spin vector more 
rigorous and gauge independent. There are some inter- 
esting questions deserving further attention. Is it possi- 
ble to meaningfully compare the spin vector defined at 
the black hole with the vector defined at spatial infin- 
ity? What about comparing the spin vector of the final 
black hole with that for the individual black holes that we 
start with? Docs it matter that they are all calculated on 
different dynamical horizons? With regard to the latter 
question, it is possible that the dynamical horizons for 
the individual black holes are smoothly connected with 
the final dynamical horizon through the appearance of 
marginally trapped surfaces lying between the common 
outer horizon and the two inner horizons. This scenario 
is suggested by initial numerical studies (5?j], but fur- 
ther analytical and numerical work is required to confirm 
this and to fully understand the dynamics of marginally 
trapped surfaces. 

We also revisited the question of tidal locking lead- 
ing to corotation in close black hole binaries. Corotation 
implies both that the spin directions are aligned with 
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the orbital angular momentum and that the horizon fre- 
quency (the horizon frequency is the angular speed of 
locally-non-rotating observers as they pass through the 
horizon, as seen by stationary observers at infinity) is 
equal to the orbital frequency. We found in a previous 
study [34[ that the spin-up of the holes is too small (by 
two orders of magnitude) to reach corotation. In this 
paper we observed signs of alignment of the spin with 
the total angular momentum. The variations in the di- 
rection observed are a few degrees during the last two 
orbits. Further numerical and analytic (higher PN or- 
der) studies will be needed to see if this effect is strong 
enough during the slow inspiral phase to drive the binary 
toward spin-orbit alignment. 

Both the final magnitude, and the final direction, of 
a black-hole binary's remnant spin are of astrophysical 
interest: the former determines the efficiency of gravita- 
tional accretion, and the latter is reflected in the orien- 
tation of the inner accretion disk and (indirectly) in the 
launching direction of a jet. Our simulations are the first 
to follow the time dependence of the spin orientations in 
black hole mergers with initially mis-aligned spins, and 
the first to verify the spin-flip phenomenon: the sudden 
reorientation in spin axis that takes place when the bi- 
nary's orbital angular momentum is converted into spin 
in the final stages of the merger [l(| ■ In addition to influ- 
encing the gravitational wave forms, the spin evolution 
would also be reflected in any electromagnetic signature 
due to gas in orbit around the black holes. Predicting 
the latter signature is beyond the scope of the present 
paper but is a fruitful topic for further study. 
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APPENDIX A: PREVIOUS STUDIES 

In a previous paper we reported the first fully- 
nonlinear studies of highly-spinning black-hole bina- 
ries [35j], where we found that the spin can profoundly 
affect the orbital dynamics of the last pre-merger stages. 
In Ref. [35l | we studied cases where the spins were aligned 



or counter aligned with the orbital angular momentum. 
As a result of the spin-orbit coupling the merger times 
dramatically changed with respect to the non-spinning 
case. For example, for initial data corresponding to a 
quasi-circular orbit with period T ~ 125M and orbital 
frequency to = 0.05/M, the non-spinning holes would or- 
bit twice before merging into a single horizon, while the 
spinning holes aligned with the orbital angular momen- 
tum and spinning at a rate S/m = 0.75 (where S is 
the magnitude of the spin angular momentum and m is 
the mass of the black hole) would orbit three times be- 
fore merger. The anti-aligned spinning holes with specific 
spins S/m 2 = —0.75 would only complete one orbit be- 
fore the common event horizon formed. These results can 
be summarized by a linear fit to the Richardson extrap- 
olated merger times tcAH (formation time of the first 
common apparent horizon) of the most accurate runs 
with S/m 2 = 0.0, S/m 2 = 0.1, and S/m 2 = -0.757 
(see Refs. Hill) 

= (172 ± 1) + (40 ± 2) S, (Al) 

M 

where S = {S\/m\ + Si/m%)i. Note that extrapolating 
to maximally spinning holes gives a merger time (from 
orbital u = 0.05/Af) of 87 M and 255Af for anti-aligned 
and aligned spins respectively. 

Extrapolation to maximally rotating black holes 
aligned with the orbital angular momentum leads to rem- 
nant black holes having a sub-maximal specific rotation 
parameter S/m 2 < 0.95 which implies one cannot gen- 
erate extreme rotating black holes or violate the cosmic 
censorship hypothesis starting from orbiting black holes 
(see Fig. 02]). 

A quadratic fit to the remnant black hole of the merger 
of aligned or anti-aligned spinning holes produces 

{S/Mlj)\ R = 0.6879 + 0.1476 (5) - 0.00935 (S) 2 , (A2) 

while a fit to the energy radiated versus the initial indi- 
vidual spins yields 

-^f = 0.0348 + 0.01485 (S) + 0.00425 (S) 2 , (A3) 

While we expect more simulations of spinning black 
holes for other values of the individual spins and with 
even higher accuracy will give improved fits, Eqs. (|A1[) - 
(|A3|) already provide valuable information for data ana- 
lysts and for theoreticians modeling the merger of spin- 
ning black-hole binaries with post-Newtonian or 'Kludge' 
waveforms 63]. 

This differential orbital dynamics in turn also notably 
changes waveforms (see Figs. 1-3 in Ref. (35l]V 

We then explored changes in the magnitude of the spin 
due to tidal effects in binaries and the transfer of orbital 
angular momentum to spin and vice versa [34j . Those 
studies concluded that it is very unlikely that black-hole 
binaries become tidally locked in a corotating state dur- 
ing the last orbital stages. We considered two representa- 
tive cases, one starting with initially non-spinning black 
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holes and tracked the spin-up during the last two orbits 
before merger. The second simulation began with the bi- 
nary in an instantaneously corotating state at the same 
starting point and again tracked the spin-up of the indi- 
vidual holes. In both cases the spin-up was two orders of 
magnitude smaller than that needed to lock the binary 
into a corotating state. 



FIG. 21: A linear-least squares fit of the remnant Kerr 
spin parameter (left y-axis) and radiated energy (right y- 
axis) for the merger of equal-mass equal-spin binaries with 
spins pointing along (or in the opposite direction to) the or- 
bital angular momentum. The fits have the functional form 
y = co + ci(S/m 2 ) + c 2 (S/m 2 ) 2 . 
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